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ABSTRACT 

The  almost  exponential  increase  in  computer  power 
over  the  last  few  decades  has  opened  up  the  venue  for 
increasingly  more  advanced  and  accurate  computer 
simulations.  As  little  as  ten  years  ago,  quantum 
mechanical  calculations  were  restricted  to  predictions  of 
static  properties  of  systems  containing  tens  of  atoms, 
thus  limiting  first  principles  explorations  to  gas  phase 
chemical  and  physical  processes.  With  today’s 
computers,  quantum  mechanical  calculations  can  easily 
be  performed  for  solids  and  liquids,  thus  opening  up 
exploration  into  condensed  phase  physico-chemical 
processes.  In  this  work,  we  demonstrate  the  ability  of 
quantum  mechanical  approaches,  in  particular  the 
density-functional  method,  to  predict  shock  properties  of 
condensed  phase  energetic  materials. 

1.  INTRODUCTION 

With  the  emergence  of  teraflop  supercomputing,  the 
inclusion  of  atomistic  simulation  methods  in  critical 
energetic  materials  (EM)  research  and  development 
programs  supporting  Future  Force  weaponry  has  become 
a  reality.  Modeling  using  atomistic  methods  provides 
detailed  information  about  the  fundamental  chemical  and 
physical  processes  that  control  the  conversion  of  an 
energetic  material  to  products,  thus  providing  the 
designer  with  information  with  which  he  can  tailor  a 
material  for  specific  performance.  Such  predictive 
capabilities  also  allow  for  screening  of  notional  materials 
upon  conception,  allowing  for  elimination  of  poor 
candidates  before  investing  in  synthesis  and  testing.  At 
this  time,  the  most  popular  atomistic  simulation  method 
is  that  of  classical  molecular  dynamics  (CMD),  which 
involves  integrating  the  classical  equations  of  motion  to 
generate  temporal  profiles  of  atomic  positions  and 
velocities,  thus  providing  a  dynamic  description  of  the 
system.  Also,  thermodynamic  information  can  be 
obtained  by  averaging  properties  evaluated  at  each 
integration  step  over  the  duration  of  the  simulation. 

The  most  significant  limitation  of  a  CMD  simulation 
is  its  dependence  on  a  model  description  of  the 
interactions  among  all  particles  in  the  system.  Typically, 
the  models  are  empirical  potentials  fitted  using 
experimental  information  and  available  quantum 
mechanical  (QM)  results.  If  the  function  accurately 
describes  the  region  of  phase  space  to  be  sampled  during 


the  CMD  simulation,  then  the  results  will  be  quite 
reliable;  if  not,  then  the  results  cannot  be  assumed  to  be 
predictive.  It  is  virtually  impossible  to  accurately 
parameterize  empirical  potentials  to  describe  interatomic 
interactions  at  all  conditions  of  interest.  In  fact,  there  are 
no  CMD  models  available  that  accurately  predict 
chemical  and  physical  properties  of  EMs  under  all 
conditions.  Therefore,  we  are  currently  faced  with  either 
enhancing  existing  imperfect  empirical  models  of  EMs 
(by  making  them  system-specific  and  thus  not 
necessarily  transferable  to  other  systems)  or  using  first- 
principles  QM  methods  which  are,  by  their  nature, 
transferable.  Since  the  latter  is  obviously  more  desirable, 
we  have  directed  a  great  deal  of  effort  toward  integrating 
QM  predictions  of  static  and  dynamic  properties  of  EMs 
into  our  suite  of  modeling  capabilities. 

To  demonstrate  this  advanced  modeling  capability 
in  an  EM  application,  we  present  the  first  QM  calculation 
of  the  shock  Hugoniot  of  an  explosive.  The  Hugoniot 
curve,  which  characterizes  a  shocked  material,  can  be 
used  in  the  design  of  new  materials  or  weaponry.  The 
system  we  have  chosen  is  nitromethane,  a  prototypical 
energetic  material  for  which  numerous  experimental  data 
are  available,  including  that  of  the  shock  Hugoniot. 
Additionally,  thermochemical  and  CMD  predictions  of 
the  shock  Hugoniot  using  empirical  methods  are 
available  for  comparison. 

CMD  calculations  of  the  shock  Hugoniot  of 
nitromethane  have  been  performed  using  an  empirical 
interatomic  interaction  potential  developed  by  Sorescu  et 
al.  (1999);  however,  this  same  model  fails  to  predict 
properties  of  certain  EMs  under  conditions  of  extreme 
compression  (Sorescu  et  ah,  2000).  In  sharp  contrast,  the 
more  accurate  and  transferable  QM  methods,  in 
particular  density  functional  theory  (DFT)  (Hohenberg 
and  Kohn,  1964)  can  readily  explore  the  extreme 
pressure  and  temperature  regimes  associated  with 
detonation  and  is  not  system  specific.  Its  main  limitation 
is  that  of  computational  expense,  which  can  be 
prohibitive  for  large  systems.  To  present  date  the  only 
other  materials  whose  shock  Hugoniot  has  been 
computed  by  DFT  are  nitrogen  (Mattson  2003)  and 
deuterium  (Lenosky  et  ah,  2000;  Desjarlais,  2003;  Bonev 
et  ah,  2004),  neither  of  which  are  conventional 
explosives. 


1 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1 .  REPORT  DATE  2.  REPORT  TYPE 

01  NOV  2006  N/A 

3.  DATES  COVERED 

4.  TITLE  AND  SUBTITLE 

Using  Quantum  Mechanics  to  Predict  Shock  Properties  of  Explosives 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

U.S.  Army  Research  Laboratory,  Aberdeen  Proving  Ground,  MD 
21005-5066 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release,  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

See  also  ADM002075. 

14.  ABSTRACT 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF:  17.  LIMITATION  OF 

18.  NUMBER  19a.  NAME  OF 

a.  REPORT  b.  ABSTRACT  c.  THIS  PAGE  |J|J 

unclassified  unclassified  unclassified 

6 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


As  an  exploratory  study  on  the  suitability  of  DFT  to 
predict  shock  Hugoniots  for  materials  of  interest  to  the 
Army,  we  have  first  generated  the  shock  Hugoniot  of 
diamond.  This  system  was  chosen  because  of  the 
availability  of  experimental  data  and  theoretical  models, 
both  at  the  classical  and  quantum  mechanical  level.  The 
comparison  is  shown  in  Figure  1  (experimental  data  is 
taken  from  Pavloskii  (1971)  and  Kondo  et  al.  (1983)).  As 
evident  in  this  figure,  the  empirical  potential  of  Tersoff 
(1988)  demonstrates  unphysical  behavior  at  relatively 
small  degrees  of  compression,  while  DFT  concretely 
reproduces  experimental  information.  Thus,  these 
calculations  clearly  demonstrate  the  superiority  of  DFT 
for  use  in  predicting  the  shock  Flugoniot  of  materials. 
This  results  presented  in  this  study  will  provide 
motivation  for  future  investment  into  development  and 
integration  of  accurate  quantum  mechanical  methods  into 
the  suite  of  tools  for  computational  modeling  of  EMs. 

Materials  of  interest  to  Future  Force  weapons  (e.g. 
energetic  materials,  advanced  armor)  are  destined  to  be 
subjected  to  extreme  impacts  leading  to  substantial 
physical  or  chemical  changes  in  the  atomic  structure  of 
the  material.  Molecular  dynamics  (MD)  simulation,  an 
atomistic  simulation  method  that  provides  a  temporal 
description  of  the  behavior  of  a  system  subjected  to 
initiating  events,  is  a  particularly  effective  means  to 
examine  complex  changes  occurring  in  condensed  phase 
materials  after  shock,  including  reaction  and  failure 
waves,  phase  transitions,  and  defect  formation  and 
propagation.  Such  a  detailed  examination  cannot  be 
performed  experimentally;  however,  the  detailed 
information  obtained  from  these  simulations  can  easily 
be  used  to  either  design  experiments  that  will  transform 
the  material  to  a  desired  state  upon  shock,  or  to  design 
new  materials  with  tailored  shock  properties. 
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Figure  1:  Comparison  of  the  diamond  shock  Flugoniot 
between  experiment,  the  empirical  potential  of  Tersoff 
and  DFT.  The  experimental  data  is  from  Pavloskii 
(1971)  and  Kondo  et  al.  (1983). 


2.  SHOCK  HUGONIOT 

The  equation  of  state  (EOS)  reveals  the  physical 
properties  of  a  material  given  a  set  of  thermodynamic 
conditions.  It  can  be  used  to  determine  such  properties  as 
the  heat  capacity,  bulk  modulus,  and  thermal  expansion 
coefficients.  In  general,  the  EOS  is  too  complicated  to  be 
derived  using  statistical  mechanics  except  for  simple  toy 
models  like  the  ideal  gas.  In  mose  cases,  the  EOS  is  a 
complicated  analytic  model  which  is  fit  to  data  obtained 
from  either  experiment  or  theory. 

Of  particular  interest  to  detonation  physics,  is  the 
exploration  of  the  EOS  at  extreme  conditions  of  pressure 
and  temperature.  These  extreme  conditions  can  be 
experimentally  achieved  either  under  hydrostatic  or 
shock  conditions.  The  advent  of  the  diamond  anvil  cell 
(DAC),  has  allowed  experimentalists  to  reach  pressures 
(up  to  200  GPa)  previously  only  obtainable  by  shock 
methods.  Nevertheless,  this  has  not  decreased  the 
necessity  to  study  materials  under  shock  conditions  since 
these  shock  methods  (e.g.,  gas-gun,  laser-induced,  or 
flyer-plate)  are  able  to  reach  extreme  pressures  (>  200 
GPa)  and  temperatures  (>  4000  K)  which  are  not 
unattainable  without  destroying  the  diamonds  in  the 
DAC.  Furthermore,  knowledge  of  the  behavior  of  a 
material  under  shock  conditions  provides  knowledge  of 
the  states  in  the  detonation  process. 

The  thermodynamic  states  accessible  under  shock 
conditions  is  given  by  the  Flugoniot  which  is  the  locus  of 
points  in  (E,P,V)-space  satisfying  the  condition: 

e^e^^+pXk-v,)  (d 

where  E  is  the  internal  energy,  P  is  the  pressure,  V  is  the 
volume  and  the  subscripts  0  and  1  denote  the  initial  and 
shocked  state,  respectively.  This  relation  follows  from 
conservation  of  matter,  momentum,  and  energy  for  an 
isolated  system  compressed  by  a  pusher  at  a  constant 
velocity.  In  the  canonical  (NVT)  ensemble,  it  is  often 
convenient  to  define  the  Hugoniot  function  (Erpenbeck, 
1992;  Rice  etal.,  1996) 

h(T,V)  =  E(T,V)-E(T0,V0 ) 

-±[p(t,v)+p(t0,v0)Iv0-v)  (2) 

The  locus  of  states  which  satisfy  h(Th,Vh)  =  0  is  the 
principal  Hugoniot.  Since  one  does  not  know  the  points 
in  phase  space  which  satisfy  the  principle  Hugoniot 


condition  a  priori,  it  is  almost  always  necessary  to 
interpolate  using  neighboring  points 

3.  COMPUTATIONAL  METHOD 

The  principal  Hugoniot  conditions  are  determined 
using  a  method  proposed  by  Erpenbeck  (1992).  In  this 
method,  EOS  points  over  a  range  of  T,  V  and  P  are 
generated  using  molecular  dynamics.  In  this  work,  we 
will  use  canonical  molecular  dynamics  (NVT-MD)  to 
evaluate  all  EOS  points.  Each  principal  Hugoniot  point 
is  determined  through  the  following  steps.  First,  after 
selecting  a  volume  for  which  the  Hugoniot  point  is 
desired,  NVT-MD  simulations  are  performed  for  a  few 
temperatures  believed  to  correspond  to  that  for  which  the 
Hugoniot  function  Eq.  (2)  is  zero.  The  Hugoniot 
function  Eq.  (2)  is  evaluated  for  each  temperature,  and 
this  set  of  Hugoniot  points  is  fitted  to  a  linear  equation  in 
temperature.  From  this,  it  is  trivial  to  determine  the 
temperature  for  which  Eq.  (2)  is  zero  (the  Hugoniot 
temperature).  Next,  the  pressures  corresponding  to  the 
various  simulations  are  also  fitted  to  a  linear  equation  in 
temperature.  The  Hugoniot  pressure  is  determined 
simply  by  evaluating  the  fitted  equation  at  the  Hugoniot 
temperature. 

The  principal  Hugoniot  requires  accurate  calculation 
of  the  internal  energy  and  pressure  of  a  system  at  a  given 
density  and  temperature.  We  choose  to  compute  these 
thermodynamic  quantities  using  Kohn-Sham  (Kohn  and 
Sham,  1965)  density  functional  theory  (DFT;  Hohenberg 
and  Kohn,  1964)  as  implemented  in  the  VASP  (Vienna 
ab  initio  simulation  program;  Kresse  and  Hafner,  1992; 
1994;  Kresse  and  Furthmuller,  1996;  1996)  code 
developed  at  the  Technical  University  of  Vienna.  Our 
thermodynamic  quantities  are  computing  by  performing 
Born-Oppenheimer  molecular  dynamics  (BOMD)  in  the 
canonical  ensemble  using  a  Nose-Hoover  thermostat  for 
the  ions  (Nose,  1984;  1984;  Hoover,  1985).  The  core 
electrons  are  replaced  with  projector  augmented  wave 
(PAW;  Kresse  and  Joubert,  1999)  potentials.  The 
exchange-correlation  energy  used  the  generalized 
gradient  approximation  (GGA)  of  Perdew-Berke- 
Ernzerhof  (PBE;  Perdew  et  ah,  1996). 

These  DFT  calculations  have  three  major  sources  of 
convergence  errors:  energy  cutoff,  k-point  sampling  and 
finite-size.  The  VASP  code  expands  the  electronic  states 
using  a  plane-wave  basis.  The  accuracy  of  the  electronic 
states  is  controlled  by  the  plane -wave  energy  cutoff  Ec. 
For  our  calculations,  we  set  Ec  =  520  eV  which 
converges  the  total  energy  and  pressure  to  at  least  10 
meV/atom  and  2  kbar,  respectively.  The  non-local 
projectors  for  the  PAW  potentials  are  evaluated  in  real- 
space  and  optimized  by  the  method  of  Kresse 
(unpublished).  While  the  non-local  projectors  are  known 
to  be  more  accurately  evaluated  in  reciprocal-space,  the 


real-space  projectors  are  faster  to  evaluate  and  exhibit 
better  scaling  with  system  size.  In  order  to  appropriately 
balance  accuracy  with  performance,  we  optimize  the 
real-space  projectors  with  a  strict  tolerance  in  order  to 
minimize  systematic  errors  in  the  total  energies  and 
pressures. 

The  Brillouin  zone  was  only  sampled  using  the  F 
point.  As  our  system  is  a  liquid,  the  k-point  sampling  is 
not  as  critical  as  for  a  crystalline  solid.  All  the 
calculations  presented  here  are  performed  on  a  eight- 
molecule  simple-cubic  cell.  We  did  not  assess  the  finite- 
size  effect  errors  for  this  calculation,  however  it  is 
believed  that  they  should  only  be  significant  near  a  phase 
transition.  Extensive  tests  for  finite-size  effects  would 
also  be  prohibitively  computationally  expensive. 

The  other  major  source  of  systematic  error  is  the 
approximation  to  the  exchange-correlation  energy. 
Unlike  the  more  advanced  quantum  mechanical 
techniques  that  treat  the  exchange-correlation  energy 
more  accurately,  it  remains  an  uncontrolled 
approximation  in  the  density-functional  method.  Its 
limitations  are  widely  known  and  documented  (see  e.g. 
Martin  (2004)  and  references  therein).  In  particular  it  is 
known  to  fail  for  vdW-bonded  materials,  such  as  most 
organic  molecules,  water,  noble  gases  and  other  sparse 
matter.  For  energetic  materials,  it  is  known  to 
underestimate  equilibrate  lattice  constants  and  bulk 
moduli  which  are  softer  than  experiment.  Hence,  there  is 
the  possibility  an  unpredictable  error  in  the  total  energies 
and  pressures  at  low  compression  ratios.  At  high 
compression  ratios  the  errors  from  the  vdW  interaction  is 
expected  to  be  less  important,  and  this  is  what  is  found  in 
shock  Hugoniot  DFT  calculations  of  liquid  deuterium 
(Lenosky  et  ah,  2000;  Desjarlais,  2003;  Bonev  et  ah, 
2004).  The  results  presented  here  are  for  the  low 
compression  regime.  This  was  chosen  in  order  to 
compare  with  existing  experimental  data  and  to  test  DFT 
method  in  the  “worst-case”  scenario. 

Lastly,  we  also  point  out  that  the  electronic 
temperature  was  set  equal  to  the  ionic  temperature,  i.e. 
Tei  =  Tjon,  using  a  Fermi-Dirac  function  for  the 
occupation  of  the  electronic  states.  All  though  there  are 
arguments  against  this  rational  based  on  the  fact  that  the 
DFT  is  a  ground  state  theory,  some  recent  quantum 
Monte  Carlo  results  on  liquid  helium  have  suggested  that 
the  excited  state  are  approximately  correct  (Militzer, 
unpublished).  In  any  case,  given  that  the  gap  at  F  is  at 
least  2.0  eV  (at  the  reference  volume  Vo)  there  would  be 
almost  no  occupation  of  the  excited  states  below  Tci  = 
2000  K  which  is  below  the  temperature  range  we 
explored  in  this  work.  The  occupancy  of  the  excited 
states  could  occur  earlier  if  significant  band  gap  closure 
was  induced  by  either  pressure  or  dissociation.  This  was 
not  observed  in  our  calculations. 
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4.  RESULTS 

The  internal  energy  and  pressure  for  a  given  point  in 
phase  space  (T,V)  was  determined  by  a  2  ps  equilibration 
period  followed  by  statistically  averaging  their 
instantaneous  values  (E,P)  during  a  2  ps  NVT  MD 
simulation.  Our  principle  Hugoniot  was  calculated  by  a 
linear  least  squares  fit  to  h(T.Vh)  at  a  fixed  Hugoniot 
volume  Vh  for  three  to  four  different  temperatures.  It  was 
then  a  simple  matter  to  determine  Th  which  satisfies 
h(Th,Vh)  =  0.  Similarly,  the  associated  Hugoniot  pressure 
Ph  was  determined  by  a  linear  least  square  fit  of  P(T,Vh) 
for  the  same  temperature  points;  the  Hugoniot  pressures 
is  then  given  by  Ph  =  P  (Th,Vh).  To  facilitate  comparison 
with  experiment  and  CMD,  we  chose  the  same 
thermodynamic  variables  for  our  reference  state  P0,T0,Vo 

was  0.001  kbar,  0.889  cm  /gm,  and  298  K. 1  This  was 
chosen  to  coincide  with  that  used  in  the  experiments  of 
Delpuech  and  Menil  (Delpuech  and  Menil,  1983)  and 
data  from  Brennan  (unpublished)  which  was  generated 
with  CMD  using  the  Sorescu-Thompson-Rice  (SRT; 
Sorescu  et  al.,  1999)  empirical  potential. 

Our  initial  molecular  configuration  were  generated 
by  placing  eight  nitromethane  molecules  equidistant 
from  each  other  in  a  simple  cubic  cell  and  performing  a  5 
ps  NVT  MD  simulation  at  a  temperature  of  T  =  1000  K 
to  generate  a  liquid  nitromethane.  This  was  generated  at 
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the  specific  volume  of  0.795  cm  / gm.  The  initial 
molecular  configuration  at  other  specific  volumes  was 
generated  by  a  simple  rescaling  of  the  atomic  positions 
and  lattice  vectors  which  did  not  exceed  6%.  The 
aforementioned  rescaling  process  had  the  potential  of 
exciting  intramolecular  modes  in  the  nitromethane 
molecules.  For  sufficiently  long  simulation  times, 
collisions  among  nitromethane  molecules  would 
eventually  damp  at  these  modes.  However,  the  specific 
volumes  and  simulation  times  considered  here  did  not 
meet  this  criteria.  There  is  the  possibility  that  this 
introduced  an  unwanted  bias  in  our  thermodynamic 
variables,  particularly  the  pressure.  This  point  is  the 
subject  of  further  investigation. 

Figure  2  presents  our  shock  Hugoniot  nitromethane 
results.  They  are  in  reasonable  agreement  with 
experiment  and  the  SRT  data  given  the  various  sources 
of  systematic  errors  in  our  DFT  calculation  in  particular 
the  lack  of  vdW  interactions.  The  agreement  is  expected 


1  Although  the  value  of  P0  exceeds  the  accuracy  to 
which  the  pressure  can  be  reasonably  calculated  within 
any  implementation  of  DFT,  we  found  that  our  results 
were  insensitive  to  P0  for  values  up  to  several  kbar. 


to  improve  at  higher  compression  ratios  where  the 
contribution  to  the  internal  energy  and  pressure  from  the 
vdW  interaction  is  less  significant. 


Figure  2:  Comparison  of  nitromethane  Hugoniot  among 
experiment,  the  empirical  potential  of  SRT  and  DFT. 
The  experimental  data  is  from  Delpuech  and  Menil 
(1983)  The  CMD  calculation  using  the  SRT  empirical 
potential  was  from  Brennan  (unpublished). 

CONCLUSIONS 

We  have  presented  the  first  QM  calculation  of  the 
shock  Hugoniot  for  a  conventional  explosive  using  the 
first-principles  DFT  method.  These  results  are 
encouraging  in  spite  of  the  shortcomings  of  DFT  for 
these  class  of  materials  which  are  dominated  by  the  vdW 
interaction  at  low  densities.  The  agreement  with 
experiment  and  CMD  were  reasonable  and  no  significant 
differences  were  found  for  the  limited  compression  ratios 
reported  here.  Additional,  Hugoniot  data  points  are  forth 
coming. 

It  is  important  to  point  out  that  although  the 
empirical  SRT  potential  gave  similar  results  to  those 
obtained  with  DFT,  the  latter  required  no  ad-hoc 
parametrization  of  the  fundamental  interactions.  The 
former,  on  the  other  hand,  requires  substantial  effort,  a 
large  degree  of  empiricism,  and  still  the  model  does  not 
allow  for  simulation  of  complex  behavior  such  as 
chemical  reactivity.  This  makes  DFT  a  robust  and 
preferred  method  for  predicting  materials  properties  over 
a  wide-range  of  thermodynamic  materials  for  many 
elements  on  the  periodic  table.  Additionally,  because  the 
DFT  method  is  quantum  mechanical  in  nature,  there  are 
material  properties  which  can  be  calculated  that  are 
completely  inaccessible  within  the  classical  methods, 
e.g.  band  gap,  dielectric  functions,  polarizability,  etc.. 

The  main  limitation  of  the  DFT  approach  is  its 
significant  computational  cost  which  limits  the  length 
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and  times  scales  that  can  be  simulated  using  classical 
methods.  These  are  typically  two  to  three  orders  of 
magnitude  in  time  and  system  size.  The  restriction  in 
length  and  time  scales  is  an  area  that  awaits  significant 
impact  from  algorithmic  development.  However, 
continuing  advances  in  algorithmic  development  at  the 
ARL  and  other  service  and  DOE  laboratories,  coupled 
with  increasingly  powerful  computational  platforms, 
provide  great  hope  that  such  computational  limitations 
will  soon  be  eased.  At  that  point,  routine  application  of 
quantum  molecular  dynamics  simulations  of  energetic 
materials  such  as  those  presented  here  will  be  realized, 
and  dependence  on  less  accurate,  empirical  modeling 
will  be  minimized. 
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